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Relativistic expansion of a magnetized fluid 
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We study semi-analytical time-dependent solutions of the relativistic magnetohydrodynamic (MHD) equations for the fields and the fluid 
emerging from a spherical source. We assume uniform expansion of the field and the fluid and a polytropic relation between the density 
and the pressure of the fluid. The expansion velocity is small near the base but approaches the speed of light at the light sphere where 
the flux terminates. We find self-consistent solutions for the density and the magnetic flux. The details of the solution depend on the 
ratio of the toroidal and the poloidal magnetic field, the ratio of the energy carried by the fluid and the electromagnetic field and the 
maximum velocity it reaches. 
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1 Introduction 



Prendergast (2005) presented a study of time dependent, relativistic, force-free, ideal MHD in the absence 



of matter by imposing a self-similar form for the solutions of the problem. In his pioneering work he assumed 
that time-dependence appears only through a dimensionless variable which contained in addition to the 
time, the speed of light and the radial distance from the centre. Then, he found a relativistic form of the 



Grad-Shafranov differential equation (Grad and Rubin 1958, Shafranov 1958, 1966). However, he did not 



take into account any pressure due to the surrounding plasma. In this paper we built up on Prendergast 's 
work by taking into account the effect of pressure while still aiming for equilibrium solutions. The presence 
of pressure leads to extra forces and the electromagnetic field is no more force-free, but now it is the net 
force due to the pressure and the electromagnetic field that has to be zero, thus the electromagnetic field 
interacts with the plasma. This interaction leads, in the most general case, to a set of non-linear partial 
differential equations. In this work we study forms of these equations permitting analytical solutions, so 
that we can have a general picture of such systems. 



Apart from Prendergast, problems of relativistic MHD have been studied by other authors. Chiueh et al. 



(1991) studied the asymptotic behaviour of steady, fully relativistic, axisymmetric, hydromagnetic winds 



and found that the flux surfaces take the form of cylinders and parabolas around the rotation axis, Li et al. 



( 1992 ) studied self-similar solutions for relativistic winds driven by rotating magnetic fields. Contopoulos 



(1994 



1995) studied the full ideal MHD problem for steady state cold outflows and he found that the 



form of the solution depends on the amount and the distribution of the electric current, he also presented 



self-similar solutions for the same problem. Fendt (1997) and Fendt and Greiner (2001) studied force- 



free magnetospheres near rotating black holes and found strong evidence for a hollow jet structure and 



applied these solutions to galactic superluminar sources. De Villiers et al. (2005) simulated the problem of 



magnetic accretion in a rotating black hole taking into account the Kerr metric. Heyvaerts and Norman 
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(2003) found stationary solutions for axisymmetric, polytropic, unconfined, ideal MHD wind using the 
WKB method. The general motivating for these studies is relativistic outflows in the form of jets and 
winds related either to AGN or to stellar mass black holes. They are mainly numerical and as such they 
have limitations on the parameters chosen and also on the trial functions employed to solve the systems 
of the partial differential equations. Our treatment leads to an analytical solution, where the system of 
partial differential equations simplifies by the use of self-similarity and finally we only solve an ordinary 
differential equation numerically. Analytical solutions allow an easier study of the parameter space and 
a better insight on the physical behaviour, however, they require simplifications and special boundary 
conditions. 

A way to introduce time dependence in the problem is by imposing a temporally self-similar solution. 
This type of self-similarity leads to solutions which are functions of a new variable x = r x v(r,t), which 
is a product of a power of the spatial coordinate and a combination of time and the spatial coordinate. 
This method has been used in similar problems, of which the best known is the blast wave solution by 



Sedov ( 


1946 


) and 


Taylor 


(1950 



dimensional analysis of the physical quantities involved in the system. Then the equations are solved and 
provide the details of the explosion. Examples of this type of self-similarity in force-free relativistic MHD 



can be found in Prendergast (2005), Gourgouliatos and Lynden-Bell (2008), Gourgouliatos (2009). 



The other form of self-similarity we are going to use is related to the separation of variables. In problems 
depending on two spatial variables, one can seek solutions which are products of a function of the angular 
coordinate and a function which depends on the distance from the origin. If a suitable form is imposed 
for the angular function and then the equation is solved numerically for the other function we find the 
meridionally self-similar solutions. In the case of radially self-similar solutions a suitable form is imposed 
for the radial function, and then the equation is solved for the angular part of the problem. Examples 



of MHD problems solved by this form of self-similarity can be found in 


Blandford and Payne 


Lynden-Bell and Boily 


(1994 


), 


Sauty and Tsinganos 


(1994 


), Vlahakis and Tsinganos 


(1998 


)• 



In the problem we are solving in this paper, we use the self-similarity technique in two steps. In the 
initial formalism of the problem we demand that the time evolution of the fields will only appear through 
the dimensionless combination of v = r/(ct). Then, we rewrite the system of the partial differential 
equations using this new variable. We observe that the system separates by imposing meridionally self- 
similar solutions. Thus, by choosing a class of those solutions the problem reduces to the solution of 
an ordinary differential equation, which we integrate numerically. These numerical solutions depend on 
the boundary conditions and the parameters chosen. We explore the parameter space and discuss the 
significance of the parameters chosen and their implications for the nature of the system. Tsui and Serbeto] 
(2007) studied a similar problem of an expanding magnetized fluid in the non-relativistic limit while taking 



into account Newtonian gravity. 



2 Formulation of the problem 



We consider a system containing an electromagnetic field E, B as measured by an observer stationary 
relative to the centre of the system and a fluid of rest density po and pressure p. The system expands 
uniformly with scaled velocity 



ct 



(1) 



As a result of the assumed uniform expansion, the Lagrangian derivate of the velocity vanishes and thus, 
each element of the system moves with constant velocity. This is a requirement for an equilibrium expan- 
sion. Had the Lagrangian derivative not been zero then the same fluid element would have suffered some 
acceleration or deceleration and the net force would have not been zero. The second assumption is that of 
axial symmetry, thus the physical quantities do not depend on the 4> coordinate. The third assumption is 
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the ideal MHD approximation, therefore in the frame of the fluid the electric field vanishes, thus 



E = -v x B . 



(2) 



The electromagnetic field has to satisfy Maxwell's equations 



V-B = 0, 



(3) 



V x E 



1 <9B 

~c~di 



(4) 



c 



(5) 



and 



1 dE 4vr . 
VxB = -— + — j 

c at c 



(6) 



Equations ^ and ^ which contain charge and current densities, allow us to determine these densities. 
The fluid has to satisfy the baryon mass conservation which is 



+ cv • V ) (7/Oq ) + c7/>oV • v = . 



d_ 

dt 



where 7 = (1 — v 2 )- 1 ' 2 is the Lorentz fa ctor and pp is the rest mas s density. 
The momentum equation is (see, e.g., Vlahakis and Konigl 2003) 



(7) 



7Po (J^ + cv • v) (^7cv) - Vp + 



i°E + j x B 



0, 



(8) 



where the relativistic specific enthalpy (over c 2 ) for a polytrope with T = 4/3 is 



e-i+4-2,- 



(9) 



We have chosen T = 4/3 to allow self-similar solutions ( Low||l982~ ) . Finally the entropy equation is 



1— 



(10) 



The above system of equations has to be solved in order to determine the density and the fields. 

We express the magnetic field in terms of two quantities, P and T. The flux function P depends on v 
and 6, and is the magnetic flux that passes through a cap of semi-opening angle 6 and lies in distance v 
from the origin in the velocity space. The function T is related to the toroidal component of the magnetic 
field. The expression of the magnetic field that by construction satisfies ^ is 



B 



dP 



2irr 2 sin 9\d6 



.( 



9P, 



111 
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Equation @ gives the electric field 



1 / dP 
E= - 9 . (vTee + v 2 — ej. (12) 



27rr 2 sin V <9i; 

The velocity of the field lines is vp = cExB/|B 2 | and it has a component, this component is a geometrical 
effect of the expansion and the toroidal component of the magnetic field and there is no rotation of the 
central dipole. 

For axially symmetric radial flows, the component of the momentum equation Q yields that 
(j°E + j x B) • = 0, which leads to a differential equation for P and T, 

dTdP dTdP v 2 + l BP _ 

lh~dJ ~ ~de Ih ~ v(l - v 2 ) ~m~ ' (13) 
or, by multiplying (13) with (1 — v 2 )/v, 

d fl — v 2 \dP d (1 — v 2 \dP _ 

dv\ v J~dJ~d0\ v )lh = ' ( ' 

This is the Jacobian of P and [(1 — v 2 )/v]T with respect to v and 9, thus 

T = j 2 v(3(P) , (15) 

where /3(P) is an arbitrary function of P. 



3 Solution 



3.1 The entropy equation 



Equation (10) yields a relation between density and pressure p 



where Q is a function of v and 
We can use this equation to find the density as a function of the pressure po = p 3 / 4 /Q 3 / 4 . 



3.2 The baryon mass conservation equation 

Substituting the above expression of the density in Q we find that the pressure has a form that can be 
conveniently written as 

4 4 

P = Po^-, (16) 

where po is a function of v and #. 
Then, the density is given by 



,PoV /4 7 3 ^ 3 M « 
P0 = -la" • ( 17 ) 



3.3 Maxwell's equations 

By construction, the form of the magnetic field chosen satisfies ([3]). The induction equation Q is also 
satisfied for the adopted forms of the magnetic and electric fields. The other two equations have the 
current and charge densities that are not determined yet. By solving equations (Tsl) and (I6J) for j° and j 
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respectively, we express these quantities in terms of the functions P and f3. The resulting expressions for 
the charge and current densities are 

c 87r 2 r 3 dP 86 ' 1 j 

c r d/3 r 2 dP „ dP i r v 2 d 2 P ,<9P 2 „ d 2 P i i 
87r^r d sm 6 1 L dP i o0 av \ V Y ov z ov <9(cos#) - 1 J 

Next we use these results to solve the momentum equation. 



3.4 The momentum equation 

The momentum equation Q contains three terms. The first one is the inertia term fj, which is proportional 
to the derivative of the relativistic specific enthalpy £. The second term, f p , is due to the pressure gradient. 
The third term, f em is due to the electromagnetic forces. We are going to evaluate each term of this 
equation and seek analytical and semi-analytical solutions. The first term is 



2„.2P A 



(20) 



The second term, using (16), is 



4 4 

^ r Vp -4 7 V-e r . 



(21) 



From the sum f*i and f p only the first term of f p survives and the effect of inertia is cancelled by the second 
term of the pressure force. This is because we have chosen a configuration that expands uniformly and as 
such there is no acceleration on the fluid. The fi is a pseudo-force that appears because of the choice of 
the frame of reference. 
The third term is 



167r 3 r 4 sin 2 



-VP . 



(22) 



where 



T = 



v 2 d 2 P ,dP 2 „ d 2 P 

~ 2v IT +sin ^ a 

7 Z ov A ov o{cost) 



V 2+lvP dP- 



(23) 



J- = is the relativistic form of the Grad-Shafranov equation for uniform expansion in the force- free limit, 



as it was formulated by Prendergast (2005). Indeed our study in the case of negligible pressure reduces to 



this equation, whose detailed study can be found in Gourgouliatos and Lynden-Bell (2008). In the present 



paper we study cases where the fluid pressure is no longer negligible. The total force on a volume element 
due to gas pressure and electromagnetic interaction is zero. If J- ^ the electromagnetic force is nonzero; 
it is normal to the magnetic field (since VP _L B) and, depending on the sign of points towards the 
axis or in the opposite direction. 

Note that if gravity is non-negligible, a fourth term should be added on the left-hand side of the mo- 
mentum equation Q and time and space coordinates have to be modified according to the metric. This 
gravitational term is fc 



-7 pqc £Vm/j (see e.g., Mobarry and Lovelace 



1986 



Meliani et al. 2006) 



where h = (1 — rg/r) 1 ^ 2 is the redshift factor with rs = 2GM*/c 2 the Schwarzschild radius for a central 
mass M*. The appearance of the redshift factor, which is solely a function of r, makes the separation of 
variables (v , 9) impossible. For r > rg this factor can be approximated as h ~ 1 and the gravitational 
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term simplifies to 



re 



Po 



3/4 tVgmw 4 Po /4 g 3 /y x 



Q 3 / 4 r 



V4, 



1 + 



c 2 r 



(24) 



It consists of two terms, of which the first one is proportional to r~ 5 and the second one is proportional 
to r~ 6 . This combination again does not permit self-similar solutions, as all the other terms appearing in 
the momentum equation are proportional to r~ 5 . A possible way of taking partially into account gravity 
is by assuming a plasma at distances r > rg, with non-relativistic temperatures (p <C Poc 2 ) and setting 
£ ~ 1, so that in the inertia term the derivative of £ will be taken into account, but in the gravitational 
term it will be set to £ = 1. This allows separation of variables, but also adds an extra constraint on Q. 
In our solutions we decided not take into account gravity. This is not an absurd assumption, as we are 
interested in late stages of expanding systems where the fluid has reached relativistic velocities and has 
already expanded a lot so that it is not close enough to the central mass for gravity to have an important 



effect. We remark that while preparing this paper for publication, a similar study by Takahashi et at 



(2009) appeared. These authors also attempted to include gravity; in fact the gravitational term is very 



important in their approach, since its magnitude is such as to balance the other forces for a given poloidal 
magnetic field. However, they omitted a factor j£ in the gravitational term of the momentum equation, 



compare our equation (24) with the last term of their equation (2). Our approach is different: we find the 
poloidal magnetic field that corresponds to flows in which gravity is unimportant. 



We now substitute the force densities found in equations ( 20 ) - ( 22 ) in the momentum equation to find 



FVP + 16^ sin z 6>7VVp = 



(25) 



A direct consequence is that Vpo 1 1 VP, or, 



po = Po(P) 



(26) 



Equation ( 25 ) then becomes (after substituting T from 23 ) 
2v 3 dP 



, d 2 p 

dv 2 



1 



dv 



+ 



sin 



9 d 

v^m 



1 dP 

~e~de 



sin 



+ 



(1 



2^/3— + 16tt sin 9 {i _ v2)3 — = . (27) 



This is the necessary condition that the r and 9 components of the momentum equation (TsT) are both zero, 



the 4> component of the momentum equation is zero as shown by (13). All pressure and inertia effects are 



included through the last term of the previous equation, which we shall call pressure-inertia term. In the 
case po = const the fi and f p forces cancel each other, and we are back in the force-free case J- = 0. 

As explained in Appendix [Aj the only nontrivial semi-analytic solution of the previous equation corre- 
sponds to 



P = g{v)sm 2 6. 
where cq and c\ are constants. 







d/3 
dP 



c P. 



dpp 
dP 



ci 
16vr 3 



(28) 



Equations (11), (12), (16) and (17) give the expressions of the physical quantities for the self-similar 
solution 



B 



9 



cos ve r 



yg 

2vrr 2 



sin 9eg + 



9 P 
2vrr 2 P 



7 v sin tfej 



(29) 



E 



9 P 2 2 ■ 

7 v sin ueg + 



27rr 2 P 



2vrr 2 



sin 9e^ 



(30) 
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the density and pressure are given by equations (16) and (17) and we can substitute for /3 and po 



P = ± (c P 2 + Ajo) 



1/2 



Po = Poo + 



Cl . 2 

— g^sm I 



16vr 



(31) 



and a prime denotes derivative with respect to v. Here /3oo and poo are constants, and Q is a free function 
of v and 8. The /3oo should vanish so that the azimuthal component of the magnetic field remains finite 
on the axis. Thus, /3 = ±c^ 2 g sin 2 9 and Co must not be negative. On the other hand poo should be a 
non-negative constant so that the pressure does not vanish on the axis. The positive (negative) sign of C3 
corresponds to increasing (decreasing) pressure as we move away from the axis, respectively. In addition, 
the sign of c\ controls the sign of F: they are opposite since equation (25) yields F + c\ sin 2 #7 4 u 4 = 0. For 



ci > the electromagnetic force points away from the axis and the sum of inertia and pressure forces points 
toward the axis; for ci < the opposite. Note that the freedom of an additive constant in the expression of 
the pressure means that the physical behaviour of the system does not change if we assume a background 
pressure which is proportional to 7 4 /(ci) 4 . The reason is that this part of the pressure introduces additional 



terms in the inertia and pressure gradient forces which cancel each other, see equations (20) and (21) 



Then by substituting (28) into ( |27[ ), the latter becomes the following ordinary differential equation for 
9{v) 



2 // 

v g 



3„' 



2v A g 
1 



2.9 



1 



+ 



cov 2 g 



+ 



(1-^2)2 (l_ w 2)3 



0, 



(32) 



where the first three terms depend on the poloidal magnetic flux, the fourth is related to the toroidal 
magnetic field, and the fifth comes from the sum of inertia and pressure terms of the momentum equation. 



Equation ( 32 ) can be solved numerically for various values of the parameters cq and ci . These parameters 



are related to the relative importance of the toroidal field and the fluid pressure and inertia compared with 



the poloidal field, see equations (28). We present the results of the numerical solution in the next section. 



4 Results 



We have solved numerically equation (32), which is a second order ordinary differential equation. Our 
motivation is to describe a physical system where some magnetic flux emerges from a surface v = vo and 
expands uniformly within a sphere in the velocity space extending to t> max ~ 1- We normalize g(v) to the 
dimensionless function g = g/g(vo), therefore the first boundary condition is g(vo) = 1. This choice of 
normalization affects ci which is now substituted by the dimensionless ci = c\/g(vo). The other boundary 
condition comes from the fact that the flux does not go further than v max , thus it is g(v max ) = 0. Subject 
to these boundary conditions we are going to solve the equation for various combinations of the parameters 
appearing. However, before moving to these numerical solutions we investigate its asymptotic behaviour 
for v -C 1. 



4.1 Asymptotic solutions for v 1 

When we focus to the non-relativistic limit v <C 1, the differential equation simplifies. The factor 7 = 
(1 — w 2 ) -1 / 2 is close to unity, thus the equation initially reduces to 

v 2 g" - 2v 3 g + (c v 2 -2)g + ciw 4 = 0. (33) 

The Frobenius expansion at small v with the first term g oc v gives indicial equation (F + 1)(F — 2) = 
for F < 4P]The solution with F = 2 corresponds to cylindrical field lines, while the F = — 1 to a dipolar 



1 There is also the possibility F = 4, with the solution g — (ci/10)v 4 . However, this case corresponds to overcollimated field lines and 
for that reason is rejected. 
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magnetic field. Since we are interested in a physical system where the flux is generated by a central source 
and then expands we study the dipolar solution with F = — 1. 



4.2 Parameter study 

The relative importance of the physic al q uantities determining the behaviour of the fluid is parametrized 



1/2 

by the two constants appearing on (|32|). The ratio of /3 over P is ±c , therefore a larger value of 



cq corresponds to a stronger toroidal component of the field. The definition of c\, similarly, relates the 



plasma pressure-inertia to the magnetic field. The asymptotic behaviour of (32) demonstrates that these 
terms are not important if v is small and it is the poloidal flux emerging from the central dipole that 
dominates. As these terms are multiplied by powers of 7 the solution will depend strongly on both the 
choice of the parameters and on the choice of w m ax- The pressure-inertia term is multiplied by 7 s , thus 
at very high velocities it is this term that determines the solution, whereas in intermediate velocities it is 
the combination of the parameters chosen. Therefore the solution depends on three parameters Co, ci and 
Umax- In our parameter study we examine the problem for various combinations of the parameters, the 
cases studied are for cq = 0, 0.5, 1, for ci = 0, 0.01, 0.1 and for i> max = 0.95, 0.97, 0.99. We have plotted the 
expansion factor F = dlng/dlnt; for c\ = 0.1 and cq = 1 (figure [TJ. The results are similar for cq = 0,0.5. 
The value of F demonstrates the behaviour of the field lines; F = — 1 corresponds to a dipole, and F = 2 
to a cylindrical field, F = with F' > gives the x-type points and F = with F' < gives the focal 
points of the field lines. 

We have chosen relatively small values for cq and c\ so that u max can be close to unity. Had we chosen 
larger values for these parameters, the flux would have been zero before reaching relativistic velocities. 
If we continue the integration further than i> max , even for small values of cq and ci, then g will oscillate 



around zero, these oscillations correspond to closed loops causally disconnected from the base (Tsui and 



S erbeto||2007| ). 

We also study the problem for negative values of c\ , (figure [2]) this corresponds to systems where the 



pressure decreases as we move away from the axis. We find from (31) that it is essential to include a 
background pressure poo> so that none of the gas pressure or the density are negative. 

The numerical solution for any combination of parameters verifies the fact that the field behaves like a 
dipole near the origin. Then, when it approaches the upper limit i> max the field deviates from the dipolar 
structure. When i> max comes closer to unity or c\ is relatively large and positive there is more flux generated 
which forms closed loops. These loops are contained within a separatrix surface corresponding to the root 
of the expansion factor F with F' > 0; and have focal points which form a circle on the equatorial plane 
whose radius is determined by the position of the root of F with F' < 0. The physical explanation for the 
formation of these loops is that there is more electromagnetic pressure needed to force the flux to reach a 
higher velocity. A higher velocity leads to a greater Lorentz factor multiplying the rest mass of the plasma. 
Thus a small increase in the f max leads to a dramatic increase in the pressure-inertia term and since the 
flux emerging from the spherical surface is limited, more flux has to be generated somehow in order to 
balance the forces. This extra flux appears in the form of these closed loops. The pressure-inertia term has 
also an effect on the collimation of the field: as it becomes more important the field lines become parallel 
to the axis, leading to a collimated magnetic field. The collimation and the closed loops are evident in 
figure [3] where the poloidal magnetic field lines (sections of surfaces of constant flux with a meridional 
plane) are plotted. There is also a very strong field near the 'Umax* In this area the field has a very weak 
radial component whereas its 9 component is large marking the turn over of the field lines, as they are 
not permitted to exceed t; max . By comparison to the inertia- free case we have found that the addition of 
inertia has a more important effect, as we expected, because the pressure-inertia forces, parametrized by 



ci are multiplied in (32) by a factor of 7 , thus they are very sensitive to v max . 

When ci < the system behaves differently. The direction of the pressure-inertia force term is opposite 
to the direction of the electromagnetic force arising from VP, therefore it is their relative intensity that 
determines the fate of g. The details depend on the initial conditions and the parameters but there are 
two families of solutions. For relatively large |ci| the pressure-inertia term becomes strong compared to 
the electromagnetic terms of the momentum equation early enough then g increases fast and diverges to 
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Figure 1. The expansion factor F = dlng/dlnt) according to the numerical solution of equation \32\ for Co = 1 and ci = 0.1. The solid 
line is for t> max = 0.95, the dashed one for i) max = 0.97 and the dotted one for t) max = 0.99. All solutions converge to —1 for v small so 
they have a dipolar behaviour. When v ss 1 the expansion factor becomes very small and negative, so that the field lines close within the 
light sphere. In intermediate distances the ones for v max = 0.97 and i> max = 0.99 have a positive F, so the field lines have x-type points 
and then focal points. 



infinity at v = 1, this generates infinite pressure and flux, making this configuration unacceptable. On 
the other hand for relatively small | c\ | the electromagnetic term dominates over the pressure-inertia term 
in (32) and g becomes zero for v < 1, therefore the field is confined within a sphere and there are no 
infinite fields or pressure, figure |4j These fields are acceptable. 

In non-relavistic force-free MHD the field lines coincide with the current lines, whereas now with nonzero 
displacement current the picture is not so simple as there are also forces between the charge and the electric 
field, forces due to the gradient of pressure and inertia forces. However, in our uniformly expanding model 
the electromagnetic force is normal to the lines of constant P, since it is proportional to VP, and as 
VP || Vpo so do the pressure- inertia forces. The relative importance of the various forces appears in the 
momentum equation (32), in particular v 2 g" — [2u 3 /(l — v 2 )]g' — [2/(1 — v 2 )]g expresses the electromagnetic 
forces due to the poloidal magnetic field, [cqv 2 / (1 — v 2 ) 2 ]g is the electromagnetic forces due to the toroidal 
magnetic field, and civ 4 (l — v 2 ) 3 are the forces due to pressure and inertia. We plot these terms for the 
case of Co = 1, c\ = 0.1 and v ma _ x = 0.99 in figure [5] The results for the relative intensity of forces are 
similar for other combinations of the parameters. 



5 Physical quantities 

In this section we study the relation of the parameters appearing in the equations to physical quantities. 
The physical quantities we are interested in are the energy and the twist of the field lines. 



5.1 Energy 

The magnetised fluid contains energy in four forms electromagnetic, kinetic, rest mass and thermal energy. 
The energy equation written in conservative form is dT 00 /dt + V • (cT^Xj) = where T^ v is the energy 
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F -6 



-14 




Figure 2. The expansion factor F = dlng/dlni) according to the numerical solution of equation ( |32| for negative c\. Depending on the 
choice of the parameters the solution either diverges to infinity or becomes zero for v < 1. The parameters chosen are co = 0.5, 9(0.1) = 1 
and g'(0.1) = —9.97 which are the same for all curves. Then the equation was integrated for c± = —0.01 (solid line), 5i = —0.1 (dashed 
line) and ci = —0.5 (dotted line). It is evident that when the pressure dominates over the other forces it requires more flux to achieve 
equilibrium. 



momentum tensor, (see e.g., Vlahakis and Konigl 2003), which gives 



dt 



£l 2 Poc 2 -p + 



B 2 + E< 
8vr 



+ V 



CPoc 2 7 2 cv |-ExB 
4tt 



0. 



(34) 



The first two terms inside the time derivative of the above equation represent the energy density of the 
fluid. It consists of the rest jpoc 2 , kinetic (7 — 1) jpoc 2 , and thermal [4j 2 — l) p energy densities. The last 
term inside the time derivative represents the energy density of the electromagnetic field. Similarly, the 
terms inside the space derivative correspond to the fluid energy flux and the Poynting flux. We now apply 
the Poynting theorem 



d /B 2 + E 2 



+ V • ( — E x B 

V47T 



jE. 



(35) 



equation (34) yields 



0_ 
dt 



{^ 2 p c 2 - p) + V • (^qcVcv) = j E. 



(36) 



Thus the term j • E measures the energy transfer between the fluid and the electromagnetic field. In our 
model this term equals 



j • E = j • (B x v) = (j x B) • v = cv • f e] 



cy v e g' sin 2 9 dpo 



cicy^v^g' sin 2 1 
16-7r 3 r 5 



(37) 



The sign of this quantity shows the flow of the energy, when it is positive energy flows from the field to the 
fluid and vice versa. This sign depends only the product cig', since all the other terms are positive. The 
energy density of the electromagnetic field is plotted in figure [6j the energy density of the fluid depends 
on the choice of poo an d although it is related to the energy density of the electromagnetic field it cannot 
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Figure 3. Plot of the poloidal field lines for cq = 1.0, and combinations of ci and u ma x- The first row (a, b, c) corresponds to c\ = 0, 
the second (d, e, f) to ci = 0.01 and the third (g, h, i) to ci = 0.1. The first column (a, d, g) corresponds to i) ma x = 0.95, the second (b, 
e, h) to D m ax = 0.97 and the third (c, f, i) to ^ max = 0.99. The field lines near the origin have a dipole structure, but as they approach 
v = »max they close, whereas in an ideal dipole there would not be such a boundary. The external circle marks the light sphere where 
the expansion velocity formally reaches the speed of light. The dotted line in (h) and (i) is the separatrix which corresponds to the local 
minimum of the flux and encircles the closed lobes appearing, the lobes have a central focus which lies at the equatorial plane. 




Figure 4. Left: Plot of the poloidal field lines for co = 0.5, ci = —0.01 and g'(vo) = —9.97. The pressure-inertia force is weak and causes 
little modification to the solution compared to the previous ones. Right: Plot of the poloidal field lines for cq = 0.5, ci = —0.5 and 
g'( v o) = —9-97. The separatrix field line that corresponds to expansion factor F = is the dotted line. The field lines enclosed by the 
separatrix have the usual dipolar structure, however the ones that are not enclosed emerge from a monopole at v = 1 and 9 = —it and 
finish at a second monopole at V = 1 and 8 = n. Fields containing monopoles are unphysical and thus unacceptable. These monopoles 
appear because we are trying to construct a field that is impossible, as we request uniform expansion and the total force to be zero. 
Thus we need more magnetic flux to balance the very strong force due to the inertia-pressure term. Since we allow only limited flux to 
emerge from the central sphere the extra flux emerges from the poles of the light sphere. This is the reason g goes to infinity at v = 1, 
depicting the need for extra flux. 
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Figure 5. The terms of equation i |'>2[ ) for co = 1, ci = 0.1, g(0.1) = 1 and i) max = 0.99. The solid line corresponds to the \v 2 g" — 
[2d 3 /(1 — v 2 )]g' — [2/(1 — t) 2 0)]g| the absolute value of the force due to the poloidal field, the dashed line is force due to the toroidal field 
cov 2 /(l — y2 ) 2 9 an d the dotted line is cit> 4 /(l — v 2 ) 3 . When v is small the forces of the toroidal field balance the poloidal, however near 
the top it is the the forces due to the pressure that become important. 

be defined in an rigorous way. 
5.2 Twist 

We now have an expression for the fields everywhere in space. We can evaluate the twist of the magnetic 
field lines. The lines of force of the magnetic field are determined by the relation: 



dr 
B, 



rd6 



r sin 



B,< 



(38) 



We substitute the expressions for the magnetic field ( 29 ) in the equation for the field lines ( 38 ) . From the 
first equality we take gsin 2 # = Pj, where the constant Pi is the value of P for this particular field line. 



Then we equate the first and the third part of equation (38) 



dej) 



2 1/2 i 

7 c dr 



-2rf(i-iy</)V2- 



(39) 



In this expression (39 ) 7 and g depend on v whereas the integral is to be done in r, however we can change 
the variable of integration from r to v, by including the ct of the denominator in the differential. The 
limits of integration are v = t>o, the surface the field line emerges from, and the maximum distance v = Vi 
it reaches in velocity space. This upper limit of integration is the solution of the equation g(vi) = Pi. The 
twist for a given field line between two points in velocity space is constant with time as it only depends on 
v. The physical meaning of this result is that the field lines have already been twisted before the expansion 
starts and then they merely expand radially. By symmetry the total twist for the outgoing part (from vq 
to Vi) and the incoming part (from Vi to vq) is twice the integral from vq to -Uj. The total twist is 



1/2 



7 2 dw 

(l-P/gy/2- 



(40) 



The twist explicitly depends on cq as this determines the ratio of the toroidal to the poloidal field; in 
addition, the form of g that appears in the integral depends on all the parameters and the boundary 
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Figure 6. Plot of the energy density of the electromagnetic field for Co = 1. The first line (a, b, c) has c\ = 0, the second (d, e, f) 
c\ = 0.01 and the third (g, h, i) c\ = 1. The first column (a, d, g) reaches i> max = 0.95, the second column (b, e, h) w max = 0.97 and the 
third (c, f, i) u max = 0.99. 



conditions, so apart from Co the choice of both c\ and i> max have an effect on the value of the twist. The 
presence of Pi in the integral demonstrates that the twist depends also on the individual field line we choose 
to study. The field line that emerges from the pole corresponds to Pi = and the twist only depends on 
cq and i> max , as the dependence on g is annulled since it is multiplied by zero. 



6 Discussion 



In this paper we have studied the problem of the uniform expansion of a magnetised fluid. The configu- 
ration consists of a polytrope with T = 4/3 corresponding to a completely degenerate electron gas in the 
extreme relativistic limit and an electromagnetic field that satisfies ideal MHD. This value is the only one 
allowing separation of variables, we remark that the same type of polytrope was chosen by Low ( 1982 ) 
and Tsui and Serbeto (2007) who studied non-relativistic MHD flows. On this paper we study solutions 
that are connected to the origin, whereas Tsui and Serbeto (2007) study fields that form magnetic islands 
disconnected from the origin. Their study is more appropriate for isolated magnetic plasmoids that have 
been twisted and lost causal connection with the parent object. 

The constraint of uniform expansion is introduced through the use of the dimensionless variable v = r jet 
which characterises each element of the magnetised fluid, so that it moves with constant velocity and 
suffers no acceleration; thus, the net force is zero everywhere but not the distinct forces due to the fluid 
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pressure and inertia and the electromagnetic fields. In our model we take into account the electromagnetic 
interaction, the force due to the pressure gradient and the inertia force. These forces produce work so 
energy is exchanged between the electromagnetic field and the fluid. This combination of assumptions 
leads to separation of the partial differential equations of the problem and to semi- analytical solutions. 
We do not include gravity in this model. In the presence of gravity the equations do not separate unless 
we constrained our study in non-relativistic temperatures (p <C pc 2 ) and distances r > r^. 

This structure is powered by an increasing magnetic dipole at the origin which releases new flux, thus a 
small loop of current lies at r = and the current increases with time. A mechanism that may provide an 
increasing magnetic field is the Poynting Robertson battery, this mechanism has been proposed to operate 



in AGNs by Contopoulos et al. (2009). This solution has the basic properties of a collimated jet, however 
in nature the exact form of the jet shall depend on the interaction with the external medium and the 
observational results on the radiation mechanisms which are not investigated in this paper. 



Comparing our results to the ones by Takahashi et al. (2009) we find that they are similar for small v 



where the dipolar form of the solution dominates. However when v becomes large they differ significantly. 
This is mainly of the different approach, we solve the problem making the assumption that the magnetic 
flux emerges from the origin and reaches a maximum expansion velocity and then we solve for the flux 
function, whereas they impose a poloidal field and then they solve for the other components of the equation. 

This physical system covers the late stage of the relativistic expansion of a magnetised fluid that has 
been already accelerated and each shell expands with constant velocity. As the velocity is proportional to 
the distance from the origin the inner shells do not overtake the outer ones, and there are no collisions. In 
this model the inlet boundary (vq) is also moving, so this cannot be identified with a surface remaining 
still in real space, however vq can be chosen to be small and to move slowly compared to the rest of the 
configuration. This study may be applied to systems that occur after the explosion of an object which 
contained a strong magnetic field. 

An interesting property of this model is that for some sets of parameters g' becomes positive (leading 
to a positive expansion factor F = vg'/g) for intermediate velocities, meaning that the field is collimated 
there. At larger velocities the g' becomes again negative so that the lines are closed inside the light sphere. 

The demand for semi-analytical and separable solutions sets constraints in the range of physical config- 
urations we can describe. If we are to study an accelerating or decelerating expansion it is essential to give 
up the uniform expansion parameter v = r/ct. In appendix B we show that the force equation does not 
admit self-similar solutions for an arbitrary combination of r and t. 

In this paper we have found analytical solutions for a complicated relativistic MHD problem. We are 
aware that there is a gap to bridge between the observed radiation from a relativistically expanding magne- 
tised fluid and our idealised model. We suggest that our work can be used a stepping stone for future studies 
and the check of the validity of MHD simulations. The results of the simulations can be compared with 
our solutions for consistency. We also remark that these structures allow the exchange of energy between 
the field and the fluid. This is a cooling/heating mechanism for the plasma and collimation/decollimation 
of the magnetic field respectively. 
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Appendix A Solutions of equation (27) 



We seek separable solutions of equation (27), of the self-similar form P = P(a), where a = g(v) sin 2 9. By 
using a instead of 9 we may transform from the pair of independent variables (v , 9) to the (v , a). With 
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the following elementary relations valid for any function <3?, 

d${y,6) _ d<f>(v,a) g' 8<!>{v , a) 
dv dv g da 

d$(v,6) cosO d$(v,a) 

- 2— — ^ot- 



sin ( 



we may rewrite (27), after dividing with aP' , as 

2v 3 , 2 



v 2 g" 



1 — v 



79 



1 — V 



79 



v 2 g 



da 



f3 d/3 



16vT 



1 dp 



(1 - v 2 ) 2 aP' dP (1 - v 2 ) 3 P' dP 



4ff 2 P" 
1-v 2 ! 57 



+ 



4g 
1 -v 2 



9 



aP" 

1^ 



15 



(A.l) 



(A.2) 



where primes denote derivative with respect to v or a. Note that the division with aP' is possible because 
this expression cannot be zero (this would mean that P = const., or equivalently that the poloidal magnetic 
field vanishes). 
We first examine two trivial cases. 



The first trivial case corresponds to g = A = const. In this case (A.2) yields j3 = const, p$ = const, and 



P = C(l — cos9), corresponding to a pure poloidal monopolar, force-free magnetic field. 
A second trivial case corresponds to 



g = X 



l-v 2 ' 



where A = const. In this case (A.2) yields 

-A 



dpp 
da 



16vr 3 



P ^ + 4 (P') 2 + AXP'P" + AaP'P" 
a da 



This solution corresponds to cylindrical poloidal magnetic field in the d< 1 regime, however this solution 
contains magnetic monopoles at v = 1 and 9 = 0, 7r and it is unphysical. 



m 



Vlahakis and Tsinganos 1997) 



In the general case, equation (A.2) yields that there are constants cq, c\, C2, C3, such that (see Appendix B 

Ac 2 g 2 



2 // 
v g 



2v a , 2 

1 2 g ~ 1 2' 

1 — v z 1 — v z 



cov 2 g 



(1- V 2)2 (l_x,2)3 



+ c 3 



4g {vg 



1-V 2 



9 



(A.3) 



Equation (A.2) then becomes 
~{vg') 2 



4g 
1 -v 2 



aP" 



P' 



C.3 



v 2 g 



(1_„2)2 



(3 dp 
C "-aP'dP ] ~ 



_ , R 3 1 d Po 

(1 _ ^2)3 ib7F p, d p 

4g 2 / _P^ 
l-v 2 V 2 P> 



(A.4) 



i.e., it becomes a sum of four products of a function of v with a function of a. 

We distinguish two possibilities. The first is that aP" /P' — C3 7^ 0. In that case we can divide (A.4) 
by aP"/P' — C3, and find that there are constants C4, C5 and cq such that (vg') 2 /g — 4<?/(l — v 2 ) = 
c±v 2 g/(l — v 2 ) 2 + c§v / (1 — v 2 ) 3 + AcQg 2 /(l — v 2 ). However, it is highly unlikely that the last equation has 
solutions that satisfy also (A.3), except for the trivial cases g = A and g = \v 2 /{l — v 2 ) considered above. 
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So we are left with the second possibility, which is to have aP" /P' — C3 = 0, with solution^] P 1 = a c 
Equation (A.4) then becomes 



, 2\ 2 



a 



1-V 



9- 



p d/3 \ / 16vr 3 dp 

C ° " G^dP ) + V 1 ~ ~^dP 



(A.5) 



We again distinguish the following two possibilities. The first corresponds to the case where C2 or C3 
is nonzero. In this case, by dividing (A.5) by 4c3/a — 4c2 we find that [g(l — v 2 )/v 2 ] 2 equals a linear 
combination of g(l — v 2 )/v 2 and a constant. This means that g(l — v 2 )/v 2 = const., a case that we already 
considered above. The second possibility is to have 02 = 03 = 0, in which case P = a. From equation (A.5) 
we find pdp/dP = cqo and 167r 3 dpo/dP = ci, while (A. 3) gives (32). 



In this appendix we have chosen that the angular part of the solution is proportional to sin 2 9, this cor- 
responds to a dipole field. The differential operator s'm9d/dd {^gd/dO^ admits in general eigenfunctions 
in the form of s'm9dPi(cos9)/d9, where Pi is the Legendre Polynomial of order I. The case studied above 
is for I = 1. Assuming a linear form on f3(P)d(3/dP = cqP and 167r 3 dpo/d-P = c± we find that we are 



constrained only to the dipole solution. This is because (27) reduces to the following form: 



X(v) sin9 



discos 1 
d9~ 



+ ci sin z 



(1-v 



2\3 



0. 



(A.6) 



As the pressure-inertia term is multiplied by sin 2 9 the only acceptable solution is this of a dipole for 
I = 1. In the absence of pressure and inertia there are acceptable solutions of higher order multipoles, see 



appendix B of Gourgouliatos and Lynden-Bell ( 2008 ) for more details 



Appendix B Solutions for arbitrary combination of r and t 



In this appendix we study whether it is possible to have separation of variables for a dimensionless pa- 
rameter other than v = r/(ct). Let us assume that we are looking for self-similar solutions of the form 
v = r/R(t) where R(t) is an arbitrary function of t which has dimensions of length. Let us assume a 
magnetic field of the form of equation (11). We shall follow step by step the process described in section 
2, the electric field that satisfies the induction equation is E = —[R/(cR)]e r x B. The gas pressure and 
inertia forces contribute only to the r and the 9 components of the momentum equation, whereas it is only 
the electromagnetic forces that contribute to the <j) component of the momentum equation. We equate the 
4> component to zero, (fE + j x B) • = 0. Unlike equation (13) the equation we take for it is more 
complicated 



(v s R 



dTdP 



v 2 R-l)T 



dP 
d0 



0, 



(B.l) 



where R = dR/dt. On division by (v^R 2 



the first and the second term of the above equation only 



depend on v and 9, therefore we expect the third term to have no dependence on t. This clearly happens 
when R is linear function of t and this is the condition we imposed in the first section of the paper, and 



indeed when substitute that in equation (B.l) we return to equation (13) 
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